package profile;
import static java.lang.Math.abs;
import static java.lang.Math.exp;
import static java.lang.Math.pow;
import static java.lang.Math.sqrt;
import ij.IJ;

import java.util.Arrays;

import util.QualityControlHelper;

/**
 * This class extends Profile because it has extra information needed
 * for the smoothing algorithm (SRAD)
 * 
 */
public class SradProfile extends Profile{
    
    protected double[] icov = null;
    private double[] sortedICOV = null;
    private double[] icovDeviations = null;
   
    private double profileVariance;
    private double icovMedian;
    private double icovMAD;
    
    private boolean varianceIsSet = false;
    
    private SradProcessor sradProcessor = null;
    
    public SradProfile(double[] arr) {
        super(arr);
        icov = new double[arr.length];
        
        //TODO: remove this?
        getProfileMedian();

        //have to call updateICOV functions here so we can get initial epsilon value
        updateICOV();
        updateICOVMedian();
        updateICOVMAD();
        sradProcessor = new SradProcessor();
    }
    
    public SradProfile(Profile profile) {
        this(profile.arr);
    }
    

    
    /**
     * This method returns the variance of the profile.  Normally variance is calculated
     * as Var(x) = E[(X - mean)^2] however in this case we are calculating it as
     * Var(x) = E[(X - median)^2]
     * 
     * @return variance of the profile
     */
    public double getProfileVariance() {
        if(!medianIsSet) {
            setProfileMedian();
        }
        if (!varianceIsSet) {
            setProfileVariance();
        }
        
        return profileVariance;
    }
    
    private void setProfileVariance() {
        varianceIsSet = true;
        
        setMean();
        
        int sum = 0;
        for (int i = 0; i < arr.length; i++) {
            sum += (arr[i] - profileMedian) * (arr[i] - profileMedian);
        }
        
        profileVariance = sum / arr.length;
    }
    
    private void setMean() {
        int sum = 0;
        for(int i = 0; i < arr.length; i++) {
            sum += arr[i]; 
        }
        
        profileMean = (sum / arr.length);
    }
    
    public double getProfileSTD() {
        if(!varianceIsSet) {
            setProfileVariance();
        }
        
        return Math.sqrt(profileVariance);
        
    }
    
    public double getICOVMax() {
        double max = 0;
        for(int i = 0; i < icov.length; i++) {
            if(icov[i] > max)
                max = icov[i];
        }
        return max;
    }
    
    public double getICOVMin() {
        double min = Double.MAX_VALUE;
        for(int i = 0; i < icov.length; i++) {
            if(icov[i] < min)
                min = icov[i];
        }
        return min;
    }
    

    /**
     * Initializes the ICOV array
     * 
     * @param profile
     * @return
     */
    private void updateICOV() {
        for (int i = cropInfo.left; i < cropInfo.right; i++) {
                double gradMag = QualityControlHelper.gradientMagnitude(arr, i);
                double laplacian = QualityControlHelper.laplacian(arr, i);
                
                //if arr[i] is 0, we can't divide by zero so set icov[i] to 
                //zero and continue to the next value
                if (arr[i] == 0.0) {
                    icov[i] = 0.0;
                    continue;
                }
                
                icov[i] = sqrt(abs(.5 * pow((gradMag / arr[i]), 2)
                        - (1 / 16) * pow(laplacian / arr[i], 2))
                        / (pow(1 + (1 / 4) * laplacian / arr[i], 2)));
        }
    }
    
    
    /**
     * Returns the median of the ICOV of the profile
     * 
     * @return
     */
    public double getICOVMedian() {
        if (sortedICOV == null) {
            updateICOVMedian();
        }
        
        return icovMedian;
    }
    
    private void updateICOVMedian() {
        if (sortedICOV == null) {
            sortedICOV = new double[cropInfo.width];
        }
        System.arraycopy(icov, cropInfo.left, sortedICOV, 0, cropInfo.width);
        Arrays.sort(sortedICOV);
        if (sortedICOV.length % 2 == 0) {
            icovMedian = (sortedICOV[sortedICOV.length / 2] + sortedICOV[sortedICOV.length / 2 - 1]) / 2;
        } else {
            icovMedian = sortedICOV[sortedICOV.length / 2];
        }
    }
    
    /**
     * Get's the median absolute deviation of the ICOV of the profile.  
     * 
     * @return
     */
    public double getICOVMAD() {
       //median must be calculated first, so check if sortedICOV is null
        if(sortedICOV == null) {
            getICOVMedian();
        }
        
        if(icovDeviations == null) {
            updateICOVMAD();
        }
        
        return icovMAD;
    }
    
    private void updateICOVMAD() {
        if(icovDeviations == null)
            icovDeviations = new double[cropInfo.width];
        
        for(int i = cropInfo.left; i < cropInfo.right; i++) {
            icovDeviations[i] = abs(icovMedian - sortedICOV[i]);
        }
        Arrays.sort(icovDeviations);
        icovMAD = icovDeviations[icovDeviations.length/2];
    }
    

    
    public double[] getArr() {
        return arr;
    }
    
    public double[] getICOV() {
        return icov;
    }
    
    
    public void applySrad() {
        sradProcessor.applySrad();
    }
    
    public void applySrad(double epsilon, double timestep, double sensitivity) {
        sradProcessor.epsilon = epsilon;
        sradProcessor.timestep = timestep;
        sradProcessor.sensitivity = sensitivity;
        sradProcessor.applySrad();
    }
    
    public double getSradEpsilon() {
        return sradProcessor.epsilon;
    }
    
    public double getSradTimestep() {
        return sradProcessor.timestep;
    }
    
    public double getSradSensitivity() {
        return sradProcessor.sensitivity;
    }
    
    public void setSradEpsilon(Double epsilon) {
        sradProcessor.epsilon = epsilon;
    }
    
    public void setSradTimestep(Double timestep) {
        sradProcessor.timestep = timestep;
    }
    
    public void setSradSensitivity(Double sensitivity) {
        sradProcessor.sensitivity = sensitivity;
    }
    
    
    
    public class SradProcessor {

        private double epsilon;
        private double timestep;
        private double sensitivity;

        private double[] divergence;
        private double[] diffCoeff;
        
        /*
         * Originally made this call this(double, double, double) but I couldn't figure out
         * a way to pass it the starting epsilon value since getQnot() isn't static and
         * the call to the other constructor must be the first call of the constructor
         */
        SradProcessor() {
            this.sensitivity = 4.0;
            this.epsilon = getQnot()/8;
            this.timestep = 0.1;
            divergence = new double[cropInfo.width];
            diffCoeff = new double[cropInfo.width];
        }
        
        private void applySrad() {
            double qnot = getQnot();
            int i = 0;
            
            while(true) {
              updateICOV();
              updateICOVMedian();
              updateICOVMAD();
              qnot = getQnot();
              if (qnot < epsilon || i == 10000) {
                  if(i == 10000) {
                      IJ.error("Warning: Hit max number of iterations during SRAD smoothing");
                  }
                  break;
              }
              updateDiffCoeff(qnot);
              updateDivergence();
              updateProfile(timestep);
              i++;
            } 
        }
        
        
        private void updateProfile(double timestep) {
            for(int i = 0; i < arr.length; i++) {
                arr[i] = (arr[i] + (timestep/4) * divergence[i]);
            }
        }

        private void updateDivergence() {
            for(int i = cropInfo.left + 1; i < cropInfo.right - 1; i++) {
                divergence[i] = diffCoeff[i] * (arr[i+1] - arr[i]) + diffCoeff[i-1]*(arr[i-1] - arr[i]);
            }
        }

        private void updateDiffCoeff(double qnot) {
            for(int i = cropInfo.left + 1; i < cropInfo.right - 1; i++) {
                diffCoeff[i] = exp(-(pow(icov[i],2) - pow(qnot,2)) / (pow(qnot,2) * (1 + pow(qnot,2))));
            }
        }

        private double getQnot() {
            return icovMedian + (sensitivity * sensitivityFactor());
        }
        
        private double sensitivityFactor() {
            double[] tmpIcov = new double[cropInfo.width];
            for(int i = cropInfo.left; i < cropInfo.right; i++) {
                tmpIcov[i] = abs(icov[i] - icovMedian);
            }
            Arrays.sort(tmpIcov);
            return tmpIcov[tmpIcov.length/2];
        }
    }

    

}
